According to the annotation information of the AP2 domain (PF00847), 178 AP2/ERF genes were identified from the tea plant genome. The protein sequences of these genes were extracted and compared with 147 AP2/ERF proteins in Arabidopsis. Based on the domain characteristics and sequence similarity, the AP2/ERF family in tea plant was divided into five categories, namely DREB subfamily (52 members), ERF subfamily (88 members), AP2 subfamily (30 members), RAV subfamily (4 members) and Soloist (4 members). We named these genes in accordance with the family classification and genome location information, and recorded them in Table S1. Thereafter, the basic physicochemical properties of CsAP2/ERFs were analyzed. The amino acid lengths of CsAP2/ERFs ranged from 67 aa (CsSoloists-01) to 748 aa (CsAP2-29), and the protein molecular weight (MW) ranged from 74.54 kDa (CsSoloists-01) to 811.66 kDa (CsAP2-29). The theoretical isoelectric point (pI) varied from 4.49 (CsERF-36) to 10.24 (CsERF-70).

An unrooted phylogenetic tree was constructed in tea plant by using the conserved sequences of proteins in the CsAP2/ERF family (Fig. 1). Meanwhile, the DREB and ERF subfamilies were subdivided into six groups based on two main distribution methods proposed by Sakuma [6] and Nakano [7]. Table 1 summarizes the number of AP2/ERFs of tea plant, Arabidopsis [7], poplar [11] and grapevine [10]. Overall, the ERF subfamily has a numerical advantage in all species listed. According to the comparison of the distribution of subfamilies in four species, the distribution of CsAP2/ERFs was more similar to that of poplar and grapevine. Moreover, the A3 group (IVb) of DREB subfamily was missing in tea plant and grapevine, which is consistent with previous reports [10, 26].

A neighbor-joining phylogenetic tree of the AP2/ERF family in tea plant. The phylogenetic tree was constructed based on 178 conserved domain sequences of CsAP2/ERFs. The CsAP2/ERF family was allocated to five subfamilies (DREB subfamily with groups A1-A6, ERF subfamily with groups B1-B6, AP2 subfamily, RAV subfamily and Soloist), and covered with different colors

159 CsAP2/ERFs were unevenly distributed across 15 chromosomes, except for 19 genes on contigs (Fig. 2). Chr1 had the largest number of CsAP2/ERFs (26 genes), whereas Chr10 had the smallest (2 genes). CsDREBs located on every chromosome except Chr 8 and Chr 10. Apart from that, two CsRAVs were mapped on Chr6 and Chr10, as well as three CsSoloists were mapped on Chr 11, Chr 13 and Chr15. CsERFs were more likely to get clustered on chromosomes compared with other CsAP2/ERFs. Clustered CsERFs were easy to spot in Chr1, Chr5, Chr7, Chr14 and Chr15. Every chromosome contained more than two types of CsAP2/ERFs, while Chr 8 only has CsERFs.

Genome location of 159 CsAP2/ERFs on 15 chromosomes. The chromosomal positions of CsAP2/ERFs were mapped according to the genome of tea cultivar ‘Shuchazao’. The subfamilies were shown in different colors (CsDREBs in blue, CsERFs in green, CsAP2s in yellow, CsRAVs in pink and CsSoloists in purple)

We detected gene duplication events in the CsAP2/ERF family in tea plant (Fig. 3), and 72 pairs of whole genome duplication (WGD) or segmental duplication events and 13 pairs of tandem duplication events were found. Thus, WGD or segmental duplication was the main expansion pattern of the CsAP2/ERF family. Meanwhile, tandem duplication promoted the extension of CsDREBs and CsERFs. We calculated the Ka/Ks values for all paralogous genes to assess the selection pressures (Table S4). The results showed that the ratios of Ka/Ks between all paralogous gene pairs were less than one, indicating that purifying selection was dominant during the evolution of the CsAP2/ERFs. To further explore the potential evolutionary mechanisms of the CsAP2/ERF family, the synteny analysis of the AP2/ERF families of tea plant, Arabidopsis and poplar was conducted (Fig. 4). The synteny relationships are presented in Table S5. The results showed that tea plant has more AP2/ERF gene pairs with poplar (300 pairs) than with Arabidopsis (133 pairs). This result suggested that the AP2/ERF family of tea plant was evolutionarily similar to poplar.

The synteny analysis of CsAP2/ERFs. The value on each chromosome represents the chromosome length in Mega base (Mb). The gray lines indicate all synteny blocks in the genome of tea cultivar ‘Shuchazao’, and the gold lines denote the whole genome duplication (WGD) or segmental duplicated gene pairs of CsAP2/ERFs

The synteny analysis of AP2/ERFs between Arabidopsis, tea plant and poplar. Gray lines indicate all synteny blocks between tea plant and the other two species. The red lines indicate the orthologous AP2/ERFs

The CDS, UTR and introns were analyzed to characterize the gene structure of CsAP2/ERFs. The CsAP2s had the unique gene structures in the CsAP2/ERF family, and they tended to include several short tandem CDS regions (Figure S2). Comparatively, the other four subfamilies had fewer number of CDS, ranging from one to four in the majority (Figure S3-S6), and CSS0012420 (CsDREB) was one exception which had seven CDS. All CsRAVs and several CsDREBs had a long, complete CDS almost covering the whole genes. A total of 76 CsAP2/ERFs had the UTRs in their structure, but nine of them only had the 5’-UTR as well as 16 only had the 3’-UTR.

15 conserved motifs were predicted by the MEME to investigate the key motif in CsAP2/ERFs. The motif compositions were distinct in different subfamilies (Figure S2-S6). However, motif 1 was conserved in all CsAP2/ERFs except all CsSoloists and two CsERFs (CSS0005103 and CSS0037642). CsAP2s had two main composite patterns, one was motif 1-motif 3-motif 1, and the other was motif 7-motif 10-motif 2-motif 1. Based on these basic patterns, the composition of CsAP2s would add or replace some motifs. CsDREBs, CsERFs and CsRAVs had a similar motif composition, a series connection of motif 2-motif 4-motif 1. On this basis, more than half of CsDREBs added motif 6, and only one CsDREB (CSS0033589) added motif 8. Motif 11 was conserved in CsRAVs compared with motif 9. Different from CsDREBs and CsRAVs, the motifs of CsERFs were more diverse in groups B3 and B6. Motif 9 was found in group B6, while motifs 8, 13 and 14 were found in group B3. Besides, CsSoloists were half more covered with motif 12.

The PlantCARE database was exploited to analyze the cis-acting elements in CsAP2/ERFs. As the results showed in Table S6, the elements were classified into five categories: hormone response, plant growth and metabolic regulation, stress response, structural elements and transcription factor binding sites. The hormone responsive elements include five types: abscisic acid response (ABRE), auxin response (TGA-element and AuxRR-core), gibberellin response (TATC-box, P-box and GARE-motif), MeJA response (TGACG-motif and CGTCA-motif) and salicylic acid response (TCA-element). Plant growth and metabolic regulation elements contain MSA-like (cell cycle regulation), circadian (circadian control), HD-zip 1 (differentiation of the palisade mesophyll cells), ACE (light response), CAT-box (meristem expression), RY-element (seed-specific regulation) and so on. The third type is stress responsive elements, such as the wound responsive element (WUN-motif) and the low-temperature responsive element (LTR). The fourth type consists of structural elements, such as the protein binding site (Box III/HD-Zip 3) and promoter and enhancer regions (CAAT-box). Finally, the common transcription factor binding sites include the MYB binding site (MBS, MBSI and MRE) and the MYBHv1 binding site (CCAAT-box).

In CsAP2/ERFs, the most widely distributed cis-acting elements are the structural elements, which account for more than 70% of the total amount in the five subfamilies, and are as high as 80.34% in CsSoloists. In addition, plant growth and metabolic regulation elements account for more than 10% in each subfamily, the highest is 13.98% in CsERFs, followed by 13.20% in CsAP2s and 12.72% in CsRAVs. The distribution of the hormone responsive elements widely varied, ranging from 8.67% (CsRAVs) to 3.42% (CsSoloists). CsDREBs (6.31%) and CsERFs (6.51%) have similar numbers of hormone responsive elements, while CsSoloists (3.42%) have slightly less. Stress responsive elements accounted for 3.13% (CsDREBs), 3.14% (CsERFs), 4.53% (CsAP2s), 2.89% (CsRAVs), and 4.56% (CsSoloists) of the total, respectively. Among the five subfamilies, transcription factor binding sites are the least distributed, occupying only about 1% of the total cis-acting elements.

The expression profiles of CsAP2/ERFs were detected by RNA-seq, and the results were analyzed from T1 (November 1, 2021) to T13 (March 19, 2022) to clarify the role of CsAP2/ERFs during tea plant bud break, which were classified into four stages (S1: paradormancy, S2: endodormancy, S3: ecodormancy and S4: bud expansion and break) [36,37,38]. A total of 62 CsAP2/ERFs were selected by the FPKM values (FPKM > 5) (Table S7). These genes were hierarchically clustered according to the expression similarities and grouped into six expression modules, naming clusters A-F for further analysis (Fig. 5). Cluster F contained the largest number of CsAP2/ERFs (18 members). 15 genes belonged to cluster C, followed by clusters D, B and A, which contained 11, 8 and 6 CsAP2/ERFs severally. Besides, cluster E contained the last number of CsAP2/ERFs (4 members).

Heatmap of CsAP2/ERFs during the different stages of tea plant bud break. The 62 CsAP2/ERFs clustered into six groups based on their specific expressions during the four stages (S1-S4) of tea plant bud break (S1: T1-T2, S2: T3-T7, S3: T8-T9 and S4: T10-T13). The circular heatmap showed the correlation analysis between the environmental factors and gene expression levels (DMT: daily mean temperature, DHT: daily maximum temperature, DLT: daily minimum temperature). The graph on the right of the heatmap showed the expression patterns of the six distinct clusters. Gene expression levels were represented by standardized FPKM values. The standardization method was z-score, z = (x − µ)/σ (x: original value, z: transformed value, µ: mean and σ: standard deviation)

The analysis of the clustering results indicated that there were two main expression patterns. Clusters A and B more actively expressed in the stages close to bud break (S3 and S4), while other clusters in the early stages (S1 and S2). The expression of cluster A sharply decreased after S1, and it was almost not expressed in the whole S2. The expression of cluster B was similar to that of cluster A in this phase. The expression recovery of clusters A and B was observed in S3 and slightly declined in S4. In contrast, other clusters were inactive in both S3 and S4 periods except for cluster E with a transient recovery of expression in T9 (S3). Clusters C, E and F showed apparent expression peaks during the whole expression process compared with clusters A, B and D. The highest expression levels were evident in T3 (S2), T9 (S3) and T7 (S2). The expression peak of cluster E appeared at T9, and two obvious fluctuations occurred before this. Clusters C and F had virtually identical expression patterns, and they showed high expression levels at T3 and T7. In contrast with cluster C, the expression peak of T7 was higher than that of T3 in cluster F.

WGCNA was performed to explore the potential interaction genes of CsAP2/ERFs to further elucidate the mechanism of CsAP2/ERFs in tea plant bud break regulation (Fig. 6a). In the clustering module of WGCNA, the previously mentioned CsAP2/ERFs (picked by FPKM > 5) were mainly divided into three modules, namely, Blue, Brown and Turquoise. Meanwhile, some reported genes involved in the bud break of woody perennials appeared in these three modules. The co-expression network between CsAP2/ERFs and bud break related genes were analyzed (Fig. 6b). The top 50% genes in each network were selected according to the degree values to further analysis. Subsequently, referring to the correlation coefficients (Table S8), the expression profiles of 12 CsAP2/ERFs and nine highly related genes (|r| > 0.70) were shown in Fig. 6c.

Screening of potential interacting genes based on WGCNA. (a) Clustering dendrogram of the average network adjacency for identifying potential interacting genes. The genes in modules are marked with different colors. (b) Gene networks in blue, brown and turquoise modules. The CsAP2/ERFs in these three modules are severally colored with the module color, and the orange dots show the bud break related genes. The dot size represents the degree values. (c) The expression profiles of the selected genes from the blue, brown and turquoise modules

CSS0041210 and CSS0038945 was classified into Module Blue. Three cyclin-related genes (CSS0024392, CSS0007207 and CSS0012344) were found in this module. These genes had similar expression patterns with CSS0038945 (r > 0.80) and were negatively correlated with the expression of CSS0041210 (r ≤ −0.85). The expressions of genes listed in Module Brown were consistent (r > 0.70), and their expression peaks appeared at S4 and decreased with the development of the tea buds. The gene expression peak in Module Turquoise mainly appeared in S2, while CSS0022420 was not active in this phase. The results of the correlation analysis showed that the expression of CSS0022420 was negatively correlated with CSS0041853 (CO2), and the correlation coefficient was −0.73. Concurrently, CSS0041853 (CO2) was positively correlated with the expressions of CSS0010538, CSS0008086 and CSS0037896, and the correlation coefficients were 0.78, 0.86 and 0.95, respectively. The expression of CSS0010538 was also highly consistent with CSS0003691 (DAM) and CSS0033241(CUC1), with the correlation coefficients of 0.70 and 0.79 separately. CSS0033241(CUC1) had the highest expression correlation with CSS0008086 (r = 0.80).

The WGCNA analysis mentioned above found 12 CsAP2/ERFs which had the potential relationships with bud break related genes, and the temperature experiments were performed to further verify whether these genes were involved in bud break under temperature-controlled processes. And interestingly, nine CsAP2/ERFs, which responded to high (30 °C) or low temperature (4 °C) treatment, were discovered (Fig. 7a). The results of the expression levels indicated that they were all sensitive to low temperature, and CSS0041210 and CSS0049609 were simultaneously influenced by high temperature.

Expression profiles of CsAP2/ERFs under treatments. (a) Temperature treatment. MT: middle temperature (15 °C), LT: low temperature (4 °C) and HT: high temperature (30 °C). (b) Light treatment. The error bars exhibit the means ± SE (n = 3) gained from the three independent biological replicates. The letter represents the significance of the differences (LSD test, P